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Abstract 

We discuss the usage and applicability of deflation methods for the overlap lattice Dirac operator, 
focussing on calculating the eigenvalues using a method similar to the eigCG algorithm used 
for other Dirac operators. The overlap operator, which contains several theoretical advantages 
over other formulations of lattice Quantum Chromodynamics, is more computationally expensive 
because it requires the computation of the matrix sign function. The principle change made 
compared to deflation methods for other formulations of lattice QCD is that it is necessary for 
best performance to tune the accuracy of the matrix sign function as the computation proceeds. We 
present two possible relaxation strategies, one which provides a rigorous bound for the eigenvalues 
but seems to be too conservative in practice, and a second which is less conservative but, while its 
stability is not guaranteed, seems to work well in practice. 

We adapt the original eigCG algorithm for two of the preferred inversion algorithms for overlap 
fermions, GMRESR(relCG) and GMRESR(relSUMR). Before deflation, the rate of convergence 
of these routines in terms of iterations is similar, but, since the Shifted Unitary Minimal Residual 
(SUMR) algorithm only requires one call to the matrix sign function compared to the two calls 
required for Gonjugate Gradient (GG), SUMR is usually preferred for single inversions of the 
Dirac operator. We construct bounds for the required accuracy of the matrix sign function during 
the eigenvalue calculation. For the SUMR algorithm, we use the standard Galerkin projection 
to perform the deflation; while for the CG algorithm, we are able to use a considerably superior 
spectral pre-conditioner. The superior performance of the spectral preconditioner, and its need 
for less accurate eigenvalues, almost erodes SUMR’s advantage over GG as an inversion algorithm. 

We see factor of three gains for the inversion algorithm from the deflation on our small test 
lattices; we expect larger gains over the undeflated algorithms in realistic simulations on larger 
lattices and with smaller masses. There is, however, a significant cost in the eigenvalue calculation 
because we cannot relax the accuracy of the matrix sign function as aggressively when calculating 
the eigenvalues as we do while performing the inversions. This set-up cost is, however, more than 
compensated for the gain in the deflation if enough right hand sides are required. 

Keywords: Ghiral fermions. Lattice QGD 
PACS: ll.30.Rd, 11.15.Ha 


1. Introduction 

The approximate and spontaneously broken chiral symmetry is important in Quantum Chro¬ 
modynamics (QCD), the theory which describes the strong nuclear force, as it determines (to a 
large extent) the mass spectrum of the lightest hadrons. It also protects the effective Lagrangian 
describing the hadrons from picking up an additional mass renormalisation, which, in principle, 
will be sensitive to the ultra-violet cut-off. When simulating QCD numerically on a discrete space 
time lattice [l|, it is therefore advisable to break chiral symmetry as lightly as possible. While it 
is possible to remove the additive mass renormalisation by tuning the quark mass, this tuning is 
unlikely to be perfectly realised, creating a small error; the worse the chiral symmetry breaking, 
the harder it is to control this error. Equally, an action which violates chiral symmetry has many 
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additional operators in the effective Lagrangian, making the analysis of some observables more 
challenging. 

The difficulty is that simulating lattice QCD with an exact Ginsparg-Wilson chiral symmetry [ 2 , 
ii (the closest to the continuum chiral symmetry achievable on the lattice given the Nielsen- 
Ninomoya theorem i) is expensive. The simplest practical and known lattice Dirac operator with 
exac10 chiral symmetry is the overlap operator, mm 

aD[n] = ^{l+H + {l-n)-f5sigii{K)), ( 1 ) 

where sign(it') is the matrix sign function, /r/(l — /i) is proportional to the bare fermion mass, 
and we will call the Hermitian operator K the kernel of the matrix sign function. There is a great 
deal of freedom while choosing K - the constraints are that we require that D has the correct 
continuum limit (as the lattice spacing a goes to zero), no fermion doublers, and is local. A 
convenient choice is the Wilson Dirac operator. 


ciKxy — ys 


4 

^xy E((i ~ lp)^y,x+a(iUy,{x) + (1 + 'yfi)Sy^x-afiUj^{y)) 

A1=1 


( 2 ) 


with K = for m a real parameter in the range 0 < m < 2, and where Uy,{x) is the gauge 

connection, an SU( 3 ) matrix on every link of the lattice (we choose to use k = 0 . 19 , which has 
good locality properties for the Dirac operator). 7^ (/r = 1 ,... , 5 ) are the standard Hermitian 
form of the anti-commuting Dirac 7-matrices, which satisfy 7^ = 1 . 

It is easy to show that the massless overlap Dirac operator satisfies the Ginsparg-Wilson relation 


D[0]75(1 - 2 aD[ 0 ]) + 75^(0] = 0 , ( 3 ) 

which reduces to 75^(0] -I- D[0]75 = 0 , the equation that permits chiral symmetry, in the naive 
continuum limit (a —>■ 0 ). The overlap Dirac operator is also 75-Hermitian, D[/r]^ = 75D[/x]75, 
which guarantees that the eigenvalues of the Dirac operator are either real or come in complex 
conjugate pairs. This also means that 75D is Hermitian, and we shall call this the Hermitian Dirac 
operator. Given that 75sign(Ar) is unitary, the real eigenvalues of the massless operator are either 
at zero or 1, and the eigenvalues lie on a circle in the complex plane of radius 1/2 centred at 1/2. 
Furthermore, since the squared Hermitian Dirac operator commutes with 75, [y^jD^D] — 0 , we 
see the eigenvalues of this operator must be degenerate, and can be expressed as exact eigenvectors 
of 75 which we may denote as and ip-i, and, except when the corresponding eigenvalues are 
at 0 and ±1, the eigenvalues of 75 D occur in ± pairs which are linear combinations of ip+i and 
The eigenvectors of D are independent of the mass. Finally, since [D, D^] = 0 , the paired 
eigenvectors of D (with complex conjugate eigenvalues) are also linear combinations of ip+i and 
ip-i- As D is a normal operator (being shifted unitary), its left and right eigenvectors are the 
same, and its eigenvectors are orthogonal. 

The difficulties with using overlap fermions are due to the matrix sign function, and these are 
both algorithmic and concern the large computational cost required. The matrix sign function is 
defined in terms of a spectral decomposition, 

sign(A:) = ^ ^/i7/|sign(A), ( 4 ) 


where tpi and are the eigenvectors and eigenvalues of K. In practice, one simulates the matrix 
sign function by deflating the smallest eigenvectors of and then use some approximation to 


^The overlap operator gives an exact chiral symmetry. However, in practice, we do not render the matrix sign 
function within the overlap operator to a perfect precision, which breaks chiral symmetry. It can, however, be 
approximated to an arbitrary precision controlled by the accuracy of the approximation to the matrix sign function 
and the precision of the floating point arithmetic on the computer. The amount of explicit chiral symmetry breaking 
can therefore be controlled, and the systematic error from this approximation reduced so that it is insignificant. 
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cover the rest of the eigenvalue spectrum. Various methods have been proposed for the approxi¬ 
mation; polynomial approximations [l^ . rational approximations [m, an approach based on the 


Lanczos algorithm [l2j| and continued fractions, which are most easily expressed in terms of a five 
dimensional representation ^ 14, 15, 161. One may also use the approximation for the bulk of the 
eigenvalue spectrum without the deflation of the small eigenvectors giving a continuous function 
of K where the matrix sign function breaks down for small eigenvalues oi K 17. 18. which, 

although it loses the exact chiral symmetry, still maintains a very good chiral symmetry, and 
tends to be quicker and removes various al gori thmic difficulties associated with the discontinuity 

In this work, we shall use deflation and 
allowing us to control the accuracy of the 


in the matrix sign function [2lL l22l [23 



the optimum Zolotarev rational approximation 


matrix sign function, and therefore the breaking of chiral symmetry, with the only restriction on 
the accuracy being the numerical precision of the computer. In practical applications, since the 
kernel eigenvectors only need to be calculated once for each gauge field conhguration, the cost of 
calculating these eigenvectors (for example by using a polynomial preconditioned restarted Arnoldi 
algorithm [^ 1 is negligible. However, we need not always use the full accuracy matrix sign 
function during the computation, and, since low accuracy estimations of the matrix sign function 
are considerably faster to compute, it is advantageous to relax the accuracy for almost all calls to 
the overlap operator while still maintaining the full accuracy of the final result. 

For physical applications, we need to invert the overlap operator and, for some observables, 
to calculate its eigenvectors. There are two general strategies used to invert the overlap operator. 
The first is to invert the five dimensional representation of the matrix sign function 31| ; the second 
is to use a nested inversion 32, with an inner inversion to calculate the matrix sign function 
and an outer inversion of the overlap operator. It is not yet clear which of these approaches is 
superior. In this work, we concentrate on the nested approach. 

Two of the most efficient algorithms to invert the overlap operator are Conjugate Gradient (CG) 
and Shifted Unitary Minimal Residual (SUMR). The first of these algorithms, which is used for 
the Hermitian squared operator D, is well known and understood. The second, less well known 
option, was originally proposed by Jagels and Reichel [3113, and it is specihcally designed for 
operators such as the overlap operator. In principle, they should require roughly the same number 
of iterations to converge (although in practice without preconditioning the convergence of CG 
tends to be a little faster, though not by any significant amount). However, since SUMR requires 
one call to the matrix sign function per iteration rather than CG’s two, for a single inversion of 
the overlap operator SUMR should be the preferred method (although CG tends to be better than 
two SUMRs for inverting D^D). The relative performance of CG and SUMR is shown for one 
typical configuration in figure |2l 

Both of these methods can be improved signihcantly from the naive inversion algorithm. The 
first key idea is to use as low an accuracy as possible for the matrix sign function for the bulk 
of the calculation. This can be done firstly by relaxing the accuracy of the matrix sign function 
as the inversion progresses (3 |. and secondly by using a low accuracy inversion of the overlap 
operator as a preconditioner for a high accuracy inversion of the overlap operator (giving us three 
inversions: the inner calculation of the matrix sign function; the middle inversion - the GG or 
SUMR preconditioner, and the outer inversion of the full-accuracy overlap operator) js^l- By 
further incorporating mixed precision methods, calling a single precision matrix sign function 
when its required accuracy is sufficiently low, this approach can give an order of magnitude gain 
over the naive GG and SUMR algorithms. 

Further improvements can be made by preconditioning the GG or SUMR algorithm used in 
the middle inversion. In general, adding an additional nested inverter does not help much, but 
we have found that deflating the lowest overla p ei genvalues can gain up to an additional factor of 
five, with a higher gain at smaller quark mass [35j (our results in this work, on a small lattice and 
relatively large mass, give a factor of three improvement). Our previous approach was to calculate 
the overlap eigenvectors ■i/'i with eigenvalues Xi in advance, and then to use them to construct a 
preconditioner for the CG inversion. This approach only worked for the GG algorithm, as the 
preconditioner spoils the shifted unitary structure of the overlap operator and thus means that 
the SUMR algorithm no longer works, and there is an additional set-up cost of calculating the 
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overlap eigenvalues. However, we do not need to to calculate the overlap eigenvectors to a high 
precision for the preconditioning to be effective, and there are other ways, albeit less robust, for 
adapting deflation to SUMR. 

Of course, a better approach would be to combine the inversion and eigenvalue calculation. It 
has been known how to do this for a general operator for some time (see, for example, 3^), and the 


resulting eigCG algorithm is widely used. The purpose of this work is to report on our efforts to 
adapt the eigCG approach to overlap fermions; both to create an eigGG inverter and an eigSUMR 
inverter. The principle novelty compared to other, simpler, operators is that we need to control 
the accuracy of the matrix sign function during the inversion. The optimum relaxation strategy 
controlling the accuracy of the matrix sign function during the inversion algorithm differs from 
that of the optimum strategy for the eigenvalue calculation. Our main result is to provide bounds 
to control the accuracy of the matrix sign function as the eigenvalue calculation proceeds. Our 
second result is to confirm that these methods are superior to the undeflated inversions if enough 
right hand sides are required. Thirdly, we study the effectiveness and limitations of eigSUMR as 
an eigenvalue routine. 

Unlike simpler Dirac operators, the additional cost of eigGG or eigSUMR over the deflated 
GG or SUMR algorithm is not negligible, because we cannot relax the accuracy of the matrix 
sign function as vigorously as we would prefer during the inversion. We are instead forced to 
use a relatively high accuracy matrix sign function during the inversions used to calculate the 
eigenvectors. However, we do not need a high accuracy estimate of the eigenvectors for the 
deflation to be effective, and in general two or three inversions are sufficient to calculate the 
eigenvectors to a sufficient accuracy. Although there is some significant set-up cost, the gain from 
the deflation more than compensates for this if enough inversions are required (where ‘enough’ 
depends on the lattice volume and the quark mass). 

This article is arranged as follows. In section[2]we review the SUMR inversion algorithm, while 
in section |3] we describe the relaxation and preconditioning techniques which are currently used 
to accelerate the inversion of the overlap operator. Section |4] discusses deflation of the overlap 
inversion, both the methods we have been using to deflate the GG algorithm, and our proposal to 
apply deflation to the SUMR algorithm. In section [5] we discuss how the Krylov subspace used to 
perform the inversions can also be used to calculate eigenvalues, and we present numerical results 
in section |6l We conclude in section [T] 


2. The SUMR routine 

It has been shown in [s^ and introduced to the lattice community in [s^ that it is possible to 
use a short recurrence to build up the Arnold! vectors for a Unitary or shifted Unitary system. Both 
of these are normal matrices, so the left and right eigenvectors are the same and the eigenvectors 
are orthogonal to each other. We consider a a matrix of the form A = p + U, where U is an unitary 
N X N matrix and p is a real number (implicitly multiplied by the identity operator). The overlap 
operator can be reduced to this form with a simple scaling. Now it is clear that the Krylov subspace 
for A generated from some vector b is the same as the Krylov subspace for U starting from the 
same point. We describe the theory behind the generation of the Krylov subspace for the SUMR 
algorithm in appendix [Appendix A[ The resulting algorithm, which can be used to generate an 
orthonormal basis qt, i = 1, 2,... n which spans the Krylov subspace {b, Ab, A^b ,..., is 

given in algorithm [T] The algorithm requires a second vector, g, alongside q which also spans the 
Krylov subspace. 

It is straightforward to show by induction that in perfect arithmetic, q and q should remain 
normalised to 1. However, in practice, we have found q and q can lose their normalisation after a 
number of iterations, particularly when the matrix sign function is rendered only approximately, 
and this affects our eigenvalue routine. It is therefore necessary to monitor the normalisation 
of q, and to restart the algorithm (by normalising q, setting q = q and resetting cr, 7 and the 
other parameters needed in the SUMR inversion routine in algorithm [21) if 11911 starts growing 
significantly larger or smaller than one. 
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qo = qo = b/\\b\\ 
for j in 0,1, 2, 3,4,; do 
u = Uqj 
ij = -iqj,u) 

^3 = \/l- l7ip 

qj+1 = — (u + 

^3 

9j+i +7j'7j+i 

done. 


Algorithm 1: The modified Arnoldi algorithm to generate an orthonormal Krylov basis q for 
a unitary operator U. 


This construction of the Krylov subspace can be used to create a short-recurrence inversion 
algorithm for a shifted unitary matrix which is effectively an efficient representation of the GMRES 
algorithm, with the important improvement that there is no need to store all the q vectors 34| : 
this approach has been called SUMR in [32|. Experience has taught us that SUMR requires 


only slightly more iterations than the unpreconditioned conjugate gradient (CG) algorithm to 
converge; but it has the advantage that it only requires one call to the matrix sign function per 
iteration rather than two for GGNE, so that unpreconditioned SUMR is almost twice as fast as 
unpreconditioned GGNE for a single inversion of the overlap operator (i.e. a computation of the 
propagator). 

The full SUMR algorithm for an initial guess xo is given in algorithm [2j 

As with GG, the SUMR algorithm can be trivially be adapted to give a multi-shift solver. 
However, given that it is impossible to combine this with the various preconditioning methods 
discussed below, it is not clear that it would be beneficial to do so unless a large number of 
inversions at various distinct low masses are required. 


3. Rel^Lxation and GMRESR preconditioning for overlap fermions 

There are two principle approaches for inverting the overlap operator; firstly to use a nested 
four dimensional (4D) approach, with an inner and outer inversion [321 . [33j, l35|| , and secondly to 
express the overlap operator in terms of a five dimensional (5D) operator, and invert that 5D 
operator [H Ei, El E 1^ . To our knowledge, the optimum 4D approach and the optimum 
5D approach have not been compared. In 3l|, a 5D inversion routine was compared against a 
sub-optimum 4D method (which used relaxation, but not GMRESR or deflation), and while the 
5D inverter was shown to be superior to the relaxed 4D GG inverter by a factor of 3-4; the relaxed 
4D deflated GMRESR routine shows an even bigger gain over that routine [s^, [sll . We therefore 
still consider it an open question whether the nested 4D or 5D approach is superior. It may well 
depend on the computer architecture used. 

The previous state of the art method for the nested inversion routine was developed in 


32, 3 


|35|| . There are three steps: relaxation, GMRESR preconditioning, and deflation of the inversion 
(which should not be confused with the deflation of the matrix sign function which accelerates the 
calculation of the sign function: they are separate steps using separate sets of eigenvectors), which 
is an entirely. The key is to ensure that we do not have to evaluate the matrix sign function to a 
high accuracy during the inversion. A low accuracy approximation to the matrix sign function is 
considerably faster to evaluate than a high accuracy approximation to the matrix sign function. 
The accuracy of the matrix sign function can be measured in various ways. The accuracy of the 
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r = b- Axo-,6 = ||f||;a; = xq; (p = -;t = - 

6 p 

r 

W = 0; gold = 0;p = 0](p = 0]S = 0]q = -;q = q 

X = 0;r' = 0;r = 1;7 = l;c = l;cr = 1 
for j in 0,1, 2, 3,... 


u = Uq 
7 = 

o = a/I - | 7 p 
a = 7(5 


7 = -7 


r' = ai^ + sC/p 
f = q;(^ + 

^ 

vW+^ 


CT 


r = sa — cr 

T = —cf 
T = ST 


rj = r/r 


K = 


r'/r" 


Wo/d = w 

w = ap + K{qoid - w) 
p = p + A(goi<i - Wo/d) 
X = X + ri{q — Lo) 

5 = aS 

ip = / 5 — cip 

A = ^ 

r 

(p = s(p + 0^7^^ 


<7oid = q 
q = {'yq + u)/a 
q = aq + j^q 
done. 


Algorithm 2: The SUMR inversion routine. {q,u) = q^u represents the inner product of two 
complex vectors. 
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approximation to the sign function is defined as 


r!' = max — 1, 

hG||ft|| = l 


( 5 ) 


where s is our approximation to the matrix sign function and the maximum is over all vectors b 
normalised to 1. This is impossible to measure practically, so one usually estimates it using 


V 


/ 


h^s^b 


( 6 ) 


for one particular choice of b. The accuracy can also be controlled by tuning the approximation 
used to represent the matrix sign function. For example, the Zolotarev rational approximation for 
a real variable x, 


e{x) = 


Nz 

E 

2=1 


OJ?X 




( 7 ) 


produces an |e| which oscillates between 1 + A and 1 — A for a given range a < |a;| < b. a and b are 
fixed as soon as we have selected K: b^ is the largest eigenvalue of while is the largest of the 
small eigenvalues of used in the deflation of the matrix sign function. Nz is the order of the 
rational approximation, and the coefficients and o"^, as well as the error A, can be expressed 
as functions of Nz, a and b in terms of Jacobi elliptic functions 11, H, 371. The accuracy of the 


matrix sign function, rj, is therefore A or the error in the inversion of x^ , whichever is larger. 
A can be controlled by varying the order of Nz- In general, Nz ^ 20 allows the matrix sign 
function to be calculated to double precision accuracy. The accuracy of the matrix sign function 
can be controlled accurately by varying Nz and the accuracy of the inversion of x^ + af. 

We define the residual vector for the inversion x = A~^h as 


r = Ax — b, (8) 

and write ^ = ||r||. Here A represents either the Dirac operator for SUMR or D^D for CG. The 
question then becomes how accurate do we need ry in order that ^ is smaller than our desired 
tolerance (i.e. so that the errors from the inaccuracy of the calculation of the propagator are 
negligible compared to the other errors affecting our result)? We do not need to maintain a high 
rjQ for the whole outer iteration to invert A. We need to start with an accurate matrix sign function, 
but can gradually reduce the accuracy as the inversion progresses. The optimal relaxation strategy 
depends on the inversion routine. In a Krylov subspace method, an approximation to the residual 
is calculated at every step. If the matrix A is not exactly calculated, then this residual will diverge 
from the true residual by an increasing amount. The true residual || & — can be compared 

against the computed residual ||rfc||, leading to the inequality 

\\b-AxkW < Ikfe - (b-Axk)\\ + \\rk\\. (9) 

The goal is to control the residual gap, \\rk — {b— Axk)\\, so that it remains less than the required 
accuracy of the inversion; otherwise the true residual will stabilise at a value close to the gap. At 
the jth step of the inversion, we compute the matrix sign function to an accuracy rjj. 

For the CG algorithm, with the iteration step, 


Qj-i = = D-i - ai-Wj-i 

Pj = D- + 

lj-1 


( 10 ) 


a sensible relaxation strategy is to compute the sign function in the CG algorithm to an accuracy 
e^||6||-\/C, where C = IId||~^- This means that we want to choose r]j = eA||b||-\/C, where is 

the desired relative accuracy of the outer inversion 32, ^ . 
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32 , 


In minimal residual algorithms, including SUMR, the relaxation strategy is rjj = eA||^||/|kil| 

Is^l - We call the relaxed algorithms relCG and relSUMR respectively. 

We note that our numerical experience is that it is not possible to completely relax the accuracy 
of the matrix sign function since if it is allowed to become very inaccurate then the eigenvalue 
solver no longer converges; but instead we should ensure that the approximation to the matrix 
sign function remains better than some minimal accuracy. This minimal bound ought to be tuned 
for each ensemble; we found that an accuracy of 10“^ was effective on the ensembles we used. 

The second strategy is to use a low accuracy inversion of the overlap operator as a pre¬ 
conditioner for a high accuracy inversion of the overlap operator. Because the required accuracy 
of the matrix sign function depends on the required accuracy for the outer inversion, this means 
that the bulk of the time, which is spent in the pre-conditioner, only requires a very low accuracy 
matrix sign function. We therefore now have three inversions in the system; An outer inversion 
of A, with a middle inversion of A used as a pre-conditioner for the outer inversion, and the inner 
inversion required for the matrix sign function. The preferred routine for the outer inversion is 
GMRES, since this does not require a constant pre-conditioner, and since we only require three or 
four iteration steps for this routine, the cost of storing all the Arnold! vectors is relatively small. 
We therefore use algorithmic C and U are here arrays of vectors; c, r, u and x vectors, /3 and 

a; = 0; r = 5; C = 0; 17 = 0; z = 0 
while ||r|| > e||6||; do 

solve Au = r to relative accuracy fj 
compute c with \\Au — c|| < e||l'||||u||/||z'|| 
for j in 1, 2 ,..., z; do 
/3 = (C'[j],c) 
c = c - /3C'[j] 
u = u- l3U[j] 
done 

c = c/||c|| 
u = u/\\c\\ 

C[i -I- 1] = c 
U[i -I- 1] = zz 
a = (c, r) 

X = X + au 
r = r — ac 
i = i + l 

done (11) 

Algorithm 3: The GMRESR algorithm for overlap fermions 

a scalars, z an integer. A represents the overlap operator, e is the accuracy to which we require 
the inversion, and fj is the accuracy of the pre-conditioner, which can (and should) be tuned to 
optimise the inversion. This can be used whether we use CG or SUMR to invert the overlap 
operator, and we call the resultant algorithms GMRESR(relCG) and GMRESR(relSUMRjl. In 


^GMRESR refers to an algorithm suggested in [S^l . where a GMRES inversion was used as a pre-conditioner 
for a second GMRES algorithm as an alternative to restarting. We have replaced the preconditioner with either a 
CG or SUMR inversion. 




our numerical tests, we have neglected the additional application of the overlap operator needed 
to convert the CG algorithm into a CGNE algorithm directly comparable with SUMR. 

This approach can easily be combined with a mixed precision calculation: if the matrix sign 
function is required to a worse accuracy than ^ ~ 10 “^ (in practice, the tolerance needs 

to be a little larger than this, and to find the optimal value requires some tuning), one can use 
a single precision representation of the matrix sign function; while if it needs to be better than 
this one uses the double precision matrix sign function. This allows the bulk of the work to be 
performed in single precision, giving approximately a factor of two gain over the fully double 
precision algorithm. 

The combination of GMRESR preconditioning, relaxation and mixed precision tends to give a 
factor of ten-fifteen gain over the original unimproved algorithm. 


4. Deflation 


Eurther improvements are possible by preconditioning the relGG or relSUMR pre-conditioners. 

In [ 3 ^, it was proposed to do so by building a pre-conditioner from the eigenvectors of ff = 75 D 
(see also 3^, 40, 411 for variations on this approach). 

Suppose that we have calculated approximations to the n smallest eigenvectors of H^, ipi,ip 2 , ■ ■ ■ ,ipn, 
with eigenvalues Ai, A 2 ,..., A„. The accuracy of the eigenvectors is measured using the residual 


[Vi, Vi) 


( 12 ) 


We label the Ritz estimate of the eigenvalue as [tpi, /[tpijtpi) = where V' now represents 

a guess of the eigenvector. It is possible to construct a preconditioner for the GG algorithm 
using eigenvectors calculated with a fairly low accuracy (how low an accuracy depends on how the 
eigenvectors are calculated and the simulation parameters). We build up a spectral pre-conditioner 

= l + (13) 


used on both the left and the right of where c is a constant chosen to be somewhere in 
the bulk of the eigenvalue spectrum of H (we use c = ■^). If the eigenvectors were exact, this 
would project the smallest eigenvalues of to c, and it would be equivalent to a full 

deflation algorithm, albeit with a little more cost as we have to repeatedly apply the eigenvector 
projection. There is no need to calculate the eigenvectors accurately, and our experience is that a 
relative residual, ||ri||//Xi, of around 10 “^ for the eigenvectors is usually good enough to obtain the 
optimal convergence rate (although our results in this study seem to suggest that it does not even 
need to be as good as this). If there are a number of small eigenvalues of H separated from the 
rest of the spectrum, this simple trick may accelerate the inversion a great deal, and improvements 
in the inversion algorithm of a factor of three or four are not uncommon, with better results with 
smaller masses; of course, this has to be offset against the cost of calculating the eigenvectors 
in the first place. The disadvantage with deflation methods occur on larger lattices, where since 
the density of small eigenvalues increases, the number of eigenvectors needed to have the same 
effect also increases. This is as much of a problem for the required memory to store the additional 
eigenvectors as computational time. 

The deflated GGNE algorithm can be constructed in the standard way. For a zero initial guess, 
algorithm |4] is used to find the solution x to Hx — b = Too = 0. This is more stable and robust 
than the Galerkin deflation algorithm which is commonly used [iilillii- In particular, the 
Galerkin algorithm tends to only improve the convergence rate of the algorithm up to the point 
where the inversion residual is comparable to the accuracy of the eigenvectors, while this method 
maintains a result as good as full deflation even for very low accuracy eigenvectors. The reason 
that this algorithm is not generally used in QGD applications is because of the cost of applying the 
preconditioner, which needs to be applied on every application of the operator H^; however, for 
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ro=Po = H ^b;xQ = 0 

for j in 1, 2, 3,.. .until ||rj || < ||&||e; do 


Tj = rj_i - ajZ 

3 - 

Pj = fo- + PjPj-i 
done 

X = HH-^Xj (14) 

Algorithm 4: The preconditioned CGNE algorithm. 


overlap fermions, the additional cost of the pre-conditioner is negligible compared to that of the 
matrix sign function, so we may as well use the more robust method. This also has a disadvantage 
when using the same Krylov subspace to calculate or improve the eigenvectors, as we shall do 
below, since the Krylov subspace generated is that of the preconditioned operator rather than the 
operator we want. We found that this was not significant in practice if we use the eigenvectors 
for deflation, since only a very low accuracy of eigenvectors is required for the full effect of the 
deflation. However, it does mean that it is impractical to combine this inversion routine with a 
precise computation of the eigenvalues if they are needed for purposes other than deflation. 

However, the spectral preconditioner cannot be used for SUMR since it breaks the unitary 
structure of the operator. We must therefore find some alternative. The Galerkin algorithm uses 
some orthonormal subspace constructed from N (when N is much smaller than the size of the 
matrix A) vectors Vi which are bundled together in the matrix V to project the initial guess of 
the inversion out of the subspace. So given an initial guess xq, the new initial vector used for the 
inversion algorithm is given by 


xo = xo+ Vy^^V^{b-Axo). (15) 

This has been applied to the SUMR and other routines for overlap fermions in 0 , though without 
using the optimum relaxation and preconditioning strategies. 

For SUMR, we have used a slightly modified form of this approach. In this case, we are 
interested in the inversion of the Dirac operator D, however it is more convenient to use the 
eigenvectors of the Hermitian operator 75 D. The non-zero eigenvalues of the overlap Dirac operator 
come in complex conjugate pairs. Given that D = commutes with 75 and [D, D^\ = 0, the 
eigenvectors of D come in pairs of opposite chirality, the eigenvectors of 75 D come in pairs with 
eigenvalues of opposite sign. Each of these sets of eigenvector pairs are linear combinations of each 
other; so if an eigenvector pair of 75 D is {ip+i then the the corresponding eigenvector pair of 
D is a linear combination of these two vectors. There is no difference, then, in constructing the 
starting vector from a linear combination of the n lowest eigenvectors of 75 D rather than a linear 
combination of the eigenvectors of D. Given an orthonormal set of approximate eigenvectors ipi 
of 75 D, we first of all construct a new basis = K'j/; for unitary Y where (V'^, ) is diagonal. 

This is always possible for a Hermitian operator, since {ipi, 3 ^Dipj) would also be Hermitian and 
Y would just contain the eigenvectors of this n x n matrix. 
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We then construct the starting vector for the SUMR algorithm 


xo 

j 

where aj is chosen to minimise the norm of the residual 

N 

r' = b-Y^ ajDip'^, 

which is the same as minimising the norm of 

r = ^5b-'^ 

3 


SO 

||rf = (6,&) + ^ [|a,f - m',,b)a] 

3 


( 16 ) 


(17) 


(18) 


(19) 


where f applied to a scalar indicates complex conjugation. The residual is minimised for 


m',b) 


( 20 ) 


This procedure will obviously work for any set of vectors ip'j, not just the eigenvectors. 

If y'j were an exact eigenvector of D, this is equivalent to exact deflation; we project the low 
lying eigenvectors out of the initial residual vector. However, for approximate eigenvectors of D, 
we may still obtain a similar speed up of the inversion up to the accuracy of the eigenvectors. This 
is similar to the Galerkin projection (for an initial guess i = 0, it is equivalent to the Galerkin 
projection where A = D with a right hand side D^b), but avoids the need to invert AV for 
each right hand side. 

Unlike the spectral decomposition, this method only improves the convergence of the SUMR 
algorithm up to the accuracy of the eigenvectors. For example, if we switch off GMRESR pre¬ 
conditioning, the residual is plotted against the iteration count and number of calls to the Wilson 
operator in figurelU on a 8^ x 32 lattice, both for the original undeflated algorithm and the deflated 
algorithm. It can be seen that the rate of convergence with respect to the number of iterations 
rate of the algorithm initially improves, but then reverts to the convergence of the undeflated 
system. In terms of the iteration count, there is no real difference whether we use deflation or 
not. However, there is a small improvement in the number of calls to the kernel operator. This 
is because the initial steps of the inversion are when the most accurate matrix sign function is 
required. Fewer accurate calls to the matrix sign function are required, so the overall cost in 
the routine is reduced a little. However, the gain is still much smaller than we would usually 
desire from a deflation algorithm. This picture can be improved by improving the accuracy of the 
eigenvectors. 

This problem is averted when we use GMRESR preconditioning. We only need solve the 
preconditioner for the GMRES algorithm (the middle inversion) to a low accuracy. This means that 
while calculating the preconditioner, we remain in the regime where the projection is advantageous. 


5. Calculation of Eigenvalues 

The main challenge with deflation is the calculation of the eigenvalues. One can, of course, 
use a previous calculation of the eigenvalues, using, for example, an implicitly restarted Lanczos 
method, or a GG minimisation of the Ritz functional. An alternative is to use the Krylov subspace 
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Figure 1: The convergence in terms of the number of iterations of the inversion and the number of Wilson operator 
calls (counting a double precision Wilson operator as two calls) of the a) undefiated SUMR routine; b) deflated 
SUMR routine, with 30 overlap eigenvalues calculated with residuals below 10“^ on a 8^ x 32 dynamical overlap 
lattice. 
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generated for an inversion to calculate the eig envalues. This has previously been adapted for the 
CG algorithm for non-overlap fermions in [3g; the resulting algorithm was called eigCG. 

Suppose that we want to construct n eigenvalues V'Sj • ■ The basic idea is to construct 

a basis of vectors v from the GG residuals, {vi^ 'Vj, ■ ■ ■} and the matrix Mij = {Dvi, Dvj). After m 
iterations, one diagonalises M using some unitary operator U, so WMU = A where A is diagonal 
and sorted so that its smallest element has the smallest index. The eigenvalues of M form an 
estimate of the eigenvalues of D while the vectors v' = C/v form estimates of the eigenvectors. 
We then construct a new set of v vectors with the first n elements being the first n elements of v', 
and the remaining m — n elements being built up from the next residuals of the GG inversion, after 
re-orthogonalising v if necessary (although in practice it is quicker and more robust to initially 
calculate M with the unorthogonalised v and manipulate M to simulate the orthogonalising of 
the eigenvectors). Until the first eigenvector-restart (when we calculate the eigenvectors of M), 
the method is equivalent to the Lanczos method (given that D is Hermitian); and, as long as 
the orthogonality of the GG residuals does not break down the convergence of the eigenvectors 
continues to follow that of the un-restarted Lanczos behaviour throughout the GG inversion. When 
we complete the initial inversion, we resume the inversion against the next input vector, deflating 
those eigenvectors already calculated to the desired accuracy, with the first n elements of the new 
V taken from our best estimates of the eigenvectors. Once all the eigenvalues are calculated to the 
desired accuracy, one switches off the eigenvalue part of the routine, and continues with a simple 
deflated GG. 

There is, of course, a small cost for eigGG over a standard GG inversion, but it is not large (ne¬ 
glecting, for the moment, the effects of relaxation and preconditioning). Most of the manipulations 
of the small matrix M, including finding its eigenvectors and eigenvalues, are negligible compared 
with the computations involving the full-sized matrix (such as the Dirac operator and spinor al¬ 
gebra). One needs a few additional matrix vector products to build up M and to re-orthogonalise 
V, and a few vector manipulations to reconstruct v'. For overlap fermions, this cost (if we exclude 
relaxation) is entirely negligible compared against the cost of the matrix sign function. We are re¬ 
quired to store the vectors = {Dvi, Dv 2 , Dv^ ...} as well as v = {vi^V 2 , r’a • ■ •} during the sim¬ 
ulation. This allows us to easily measure the residual of the eigenvector, U = j^Dv — {v, ')^Dv)v^ 
without having to apply any additional matrix vector products. 

For overlap fermions, since the eigenvalues of D come in degenerate pairs, it is useful to 
separate these pairs to improve the stability of the algorithm. We do this by in place of Mij = 
{Dvi,Dvj) using Mij = {Dvi,Dvj) -I- 2 ^ 6 {vi,^ 5 Dvj), where 0 < 5 <C 1 is a tuned small constant. 
This Mij remains Hermitian and positive definite (if the eigenvalues of D are A?, then Af > 4/i^ 
so the smallest eigenvalues of D + are 4^^(1 — <5), but its eigenvectors will be the 

eigenvectors of 75 !? rather than an random mixture of the eigenvector pairs of D. The eigenvalue 
pairs of M will be separated by 4/r(5e|A|; which avoids the danger of the Lanczos algorithm missing 
one of a pair of degenerate eigenvectors (This, of course, does not help resolving degenerate zero 
eigenvalues). In a GG algorithm, one applies 'y^Dj^Dp; to allow a splitting of the eigenvalue pairs 
requires that we store ^y^Dp instead. Indeed, we have used both p and ^y^Dp as separate v vectors. 
This, of course, destroys the natural orthogonality of the residual vectors. 

Our principle interest is to create an eigSUMR algorithm for overlap fermions. The construction 
follows the same principles as the eigGG algorithm, and the crucial part of the algorithm, the 
generation of the subspace, can be summarised in algorithm [5] for the first inversion. We seek to 
find the first n eigenvalues to an accuracy and use a maximum subspace size m > n. Algorithm 
0 can be used either as a stand-alone eigenvalue routine or incorporated into a SUMR inverter. 
As long as we do not need to recalculate = Dv (which may be necessary if a low accuracy sign 
function is used) or relax the accuracy of the matrix sign function the additional overhead on top 
of the SUMR inversion is negligible. 

We have explicitly re-orthonormalised the v vectors for each calculation of the eigenvectors. 
This is necessary because we find that the orthogonalisation of the vectors in the SUMR can be 
quickly lost, especially when we use a low accuracy computation of the matrix sign function. To 
improve stability and speed, we have done so using a LDU decomposition. This allows us to avoid 
necessary additional spinor operations, which are both costly and a source for the propagation of 
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go = 90 = Vll^ll 
V = 0;v^ = 0;fc = 0 
for j in 0,1, 2, 3,4,; do 
u = Uqj] 


Vk = Qj] 


vP = 


( 1 -m) , (1 + m) 


-9i 


k = k + 1 


if (fc == m); then 

Recalculate the eigenvalues following algorithm [SI 
end if 


\-ljV 

qj+i = — (w + 7jgi) 

qj+i = ajqj+-i*qj+i 

done 

Recalculate the eigenvalues following algorithmic] 
done 


Algorithm 5: The eigSUMR eigenvalue routine, where U = 75 sign(A). 


Construct Eij = {vi,Vj) 

Perform the LDU decomposition on E, giving an upper triangular matrix U 
and a lower triangular matrix L = 

Construct My = (vf ,vf) + 

Find the unitary matrix V which diagonalises 
Vo,...n-l = (C17^v)o,,,,n-l 

k = n 

for i in 0,1 ,..., n — 1; do 

rv,i = vP - Vi{vi,j5vP) 
done 

Algorithm 6: The recalculation of the eigenvalues in the eigSUMR routine. 
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rounding errors. Our method is equivalent to the Gram-Schmidt algorithm. The matrix Eij = 
(vijVj) is decomposed as E = LDU, where L is lower triangular, D is diagonal and U is upper 
triangularis In our case, since E is Hermitian and positive definite, D is just the identity matrix 
and L = W (the LDU decomposition reduces to the Cholesky decomposition). 

After each restart of the Krylov subspace (i.e. each new right hand side), we repeated the 
diagonalisation procedure in algorithm [5] using a basis constructed from the previous three sets of 
computed eigenvectors (i.e. the eigenvectors as they were calculated at the end of the inversion, 
the eigenvectors from the inversion before that, and the eigenvectors outputed from the inversion 
before that). In principle, the largest contributions of the errors to the approximate eigenvectors 
come from the next highest eigenvectors. Each guess of the eigenvectors thus contains those 
vectors we want, plus the next highest eigenvectors plus some additional noise. By performing 
this diagonalisation, we hoped to remove the next highest eigenvectors from our present guesses. 
This gave a small performance gain in the computation of the eigenvectors at a negligible additional 
cost. 

5.1. Relaxation 

The question remains as to what is the optimal strategy for controlling the accuracy of the 
matrix sign function. The accuracy is determined by how good the approximation of = Dv is. 
In the routine above, we assumed that the matrix sign function was calculated to infinite accuracy. 
This is in practice impossible. Instead, we use an approximate matrix sign function s (with s the 
exact sign function) which leads to an approximate Dirac operator D and an approximate ii^. We 
can write. 


( 21 ) 

and our goal is to keep ||<5|| sufficiently small so that it has no significant effect on the estimate 
of the eigenvalue or the residual. Our experience is that using too low an accuracy for can 
be disastrous: the estimate of the eigenvalues of M rapidly diverges from the true eigenvalue 
spectrum, and the residual of the eigenvectors correspondingly grows worse with each restart. 
While we wish to relax the matrix sign function as much as possible, we cannot relax it too much. 

The first question is how we can measure ||i5|l; and unfortunately a direct measurement requires 
an application of the matrix sign function and is therefore expensive. It is, however, possible to 
get a quick order of magnitude estimate by using 

{v, 750 ^) = (u, 75U^) + (u, 75 ( 5 ). (22) 

(u, 75 U^) is an estimate of the eigenvalue of 75 D, and this quantity is real. Therefore, the imaginary 
part of (u, ■j 5 V^), only comes from the imaginary part of {v, 75 ( 5 ), 3(u, ■J 56 ). We can write {v, 75 ( 5 ) = 
||u||||(5|| cos^e*"^, in which case 

| 3 '(u, 75 U'°)|||u||||( 5 ||| cos 6 *sin(/)| = ||(5||| cos0sin(/)|, (23) 


given that ||u|| = 1 . 

An additional estimate of the error in is to consider the vector v' = 2—=—-—. In exact 

i-fj. 

arithmetic, this should be 75 signu, so ||D'|P should be one. The inaccuracy of can be thus be 
estimated by the deviation 

Ill'll' - 1 = (T^ - (24) 


^For the general LDU decomposition, we use the convention that Da, the elements of D satisfy \Da \ = 1 while 

La = Ul, 
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{ 6 ,v^) = 0, and given that — jiv = r"", where /r is the Ritz estimate of the eigenvalue and r" 
the residual of the eigenvector, we find that {S,v) ^ ^ 2||(5||||r"|| cos^'cosc/j'/M) where 

O' and (/)' parametrise the angle between 6 and r’'. We write 


(1 + A^ 

^ 2|m| 

a =cos9' cos(j)', (25) 


so 


j=\\ 6 f-\\S\\Pa, 


(26) 


which gives 

ll^ll = ^/3a± + 47 , (27) 

where we can cheaply measure 7 and [3 (at least to a reasonable accuracy), while a is in the range 
— 1 < a < 1. Note that we require ||i5|| > 0 and ||(5|| = when /3 = 0, which means we should 
take the positive root of the solution when 7 > 0. This means that we can establish bounds on 


max(0,-i/3+i\//32+ 47 ) < \\S\\ < + 47 ) (28) 

These results will therefore give two different estimates of ||5||. Neither of these quantities are 
ideal; the first just gives a minimum bound on <5, while the second bound depends on our estimate 
of lira'll, which can deviate significantly from the true value if ||5|| is too large. However, they each 
give us some estimate of the error without having to continually apply a high accuracy matrix sign 
function. We may thus keep track of these estimates of S, and recalculate to a high accuracy 
when one of them indicates that ||(5|| might have risen sufficiently close to es, our maximum allowed 
tolerance for S (where the precise meaning of ‘sufficiently’ has to be judged by experience). In an 
eigSUMR routine, we may also measure the accuracy of more indirectly through the difference 
between the computed and exact residuals from the inversion (the exact residual is calculated 
when we restart the pre-conditioner in an GMRESR routine). 

The residual of the estimate of the eigenvector is defined (in exact arithmetic) as 

^true = I5DV - v{v, 75 T>n) = - v{v, J5V^). (29) 

In practice in inexact arithmetic, we will have 

r'" = rtrue + 75 S - v(v, 75 ( 5 ), (30) 


which gives 


(1 - «^'^)75(5) + ((1 - t^^^^)75(5,H^™e) + (<5, (1 - ^^^'^)^). (31) 

The residual gap, g, is the difference between Hr^lp and ||r(’,.^g|p, and we want to keep this below 
e^, where e is the desired accuracy for the eigenvector. We have 

5<2|K„g|||l,5|l + ||,5f (32) 

This bound gives 

ll^ll< Ve^ + lktVeP-IKuell- (33) 
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During each update of the eigenvectors, we know that 

Vi -^{UV)ijVj 

V? ^(UV).jvf = [UV^vf + {UV),,5f, (34) 

where 5^ is either the previously calculated error from the previous eigenvector-restart if Vj is one 
of the vectors held back from the previous iteration, or due to the inaccuracy of the matrix sign 
function as we generate the new Krylov subspace. Thus, neglecting any error on (UV) due to the 
fact that they have been calculated using rather than v^, the new 6 i satisfies the bound 


ll<5,|| (35) 

i=i 

This means that to avoid having to recalculate during the iteration, we may compute the 
matrix sign function to an accuracy 


II < 


+ Iktrue.olP - \Kue,o\\ 

{m - n)kmaxn<j<m,i<n> \iUV)^j\’ 


(36) 


where k is the expected number of iterations required for the eigenvector routine to converge, or 
between high accuracy recalculations of , and ||?'"j.^g Q|| is the residual of the best converged 
eigenvector (we may use the computed residual rather than the true residual here, as the dif¬ 
ference between the two will not be significant enough to affect anything), max^ \ (UV)ij\ may 
be estimated from the previous eigenvector-restarts; the general trend in the long-term is that 
once the eigenvectors start to converge this quantity decreases as the iteration proceeds, as may 
be expected since it provides the correction to the ever-improving eigenvector. However, it can 
fluctuate significantly from one iteration to the next, and for this reason we took the average 
value from the previous ten eigenvalue calculations, which worked well in practice. Of course, this 
estimate is conservative, and in practice we relaxed it by assuming that the additional errors 6 ^ 
will not all add together coherently. It is more likely that the effects of 6 ^ will partially cancel 
each other out. We assume that the effects of these vectors on the error on will resemble a 
random walk, suggesting that need only scale the accuracy of the matrix sign function according 
to the inverse square root of the number of vectors rather than the inverse of the number of vec¬ 
tors. After the initial few calculations, even this proved to be too conservative since most values 
of \UV\ were several orders of magnitude smaller than the maximum value which only affected a 
different eigenvector each time. Therefore, we scaled the bound by the square root of the number 
of eigenvectors to simulate that only one of these eigenvectors was affected. 


+ Iktrue.oP - Iktrue.oll 
y/im- n)fcmax„/<j<m,i<„ \ {UV)ij\ 

( 1 First two diagonalisations 

( y/n Subsequent diagonalisations 


(37) 


This is the expression we have used in our simulations. While, unlike (13611 . it is not guaranteed 
to avoid large errors in the eigenvectors, it nonetheless seems to work well. When using SUMR, 
it is advisable to recalculate after the completion of the eigenvector calculation to avoid errors 
entering the calculation of the residual during deflation. 

It is also possible to reduce the accuracy of e as the algorithm proceeds; starting with a high 
value of e and gradually reducing it to the desired precision. We need to explicitly recalculate 
every k iterations to maintain accuracy, and there is no cost if we adjust e before this recalculation. 

Finally, we note that 5^ may be measured when we apply the matrix sign function to Vj. If 
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u = sv^ then 


u =sv + 75(5^ 

^-1= -Pp-^2||^ II, (38) 

where in the last step we have used the normalisation ||u|| = 1. 

When using the algorithm as a stand-alone eigenvalue solver, we may freely use the bound 
dSil), adjusting the accuracy of the matrix sign function after each eigenvector-restart. When 
incorporating this into a relaxed GMRESR(SUMR) routine, there are problems because this bound 
is not the same bound as used in the inversion. There are several possibilities: 

• Only incorporate those q vectors where the calculation of the sign function satisfies the 
bound (1571) into the Krylov subspace, at the cost of a slower calculation of the eigenvectors; 

• Relax the inversion routine as normal, including all the eigenvectors, at the cost of having 
to recalculate to a high accuracy when we construct the eigenvectors; 

• Relax the inversion routine according to the maximum of the two bounds for the accuracy 
of the matrix sign function, at the cost of having a slower inversion routine. 

In the numerical simulation presented here, we applied the third of these strategies. This 
significantly increases the cost of the initial inversions which are used to calculate the eigenvalues; 
however it means that we are able to use larger subspaces to calculate the eigenvalues so that 
they converge better; and our experience is that only three or four inversions are sufficient to have 
the smallest eigenvectors converge to the accuracy which gives the optimum rate of improvement 
for the inversions. Our experience of the first strategy on smaller test lattices showed that the 
eigenvectors did not converge, as only three or four additional vectors were used for each call to the 
GMRESR preconditioner to improve the already estimated eigenvectors: not enough to make any 
serious progress. This approach therefore offered no real gain over the undeflated GMRES(SUMR) 
routine as the eigenvectors were never of a sufficient accuracy to substantially accelerate the 
inversion. The second of these strategies is obviously inferior, as we would be calculating the 
matrix sign function twice, once to a high accuracy, per eigenvector. However, whichever strategy 
we use, there is a cost to the initial inversions used to calculate the eigenvectors. However, once we 
have calculated enough eigenvectors to a good enough accuracy, we can switch off the eigenvector 
calculation and use the optimum inversion relaxation strategy. Given a typical number of right 
hand sides, the gain from the deflation far outweighs these additional costs. It is necessary to 
recalculate after k calls to the eigenvector routine, where k can be tuned to give the best 
performance. 


6. Numerical Results 


We tested the al gori thm on 8^ x 32 dynamical overlap configurations, generated with a Luscher- 
Weisz gauge action |45l.lib. ^| at /3 = 8.15 and a quark mass of ^ = 0.03 corresponding to a pion 
mass m-jra = 0.27(4) (measured via the axial correlator, ^ 460 MeV) and a lattice spacing of 
a = 0.118(1) fm (measured from a string tension of (420 MeV)^). We used one sweep of stout 
smearing on the gauge configuration, matching the set up used to generate the ensemble. The 
simulations were run on a desktop computer at Seoul National university. We only present results 
from a single configuration with zero topological charge, since we did not observe any signihcant 
difference from conhguration to configuration; except that the improvement from the deflation was 
up to a factor of three better if the Dirac operator had an exact zero mode. We will compare the 
eigSUMR, eigGG, GMRESR(relGG) and GMRESR(relSUMR) algorithms. We set k = 25, which 
allows us to run five inversions without having to recalculate and and the desired accuracy of 
the eigenvectors to e = 10“'’’. The residual of the GMRESR preconditioner was tuned so that the 
GMRESR routine should converge in three steps after the eigenvectors were calculated and two 
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Figure 2: The convergence in terms of the number of Wilson operator calls (counting a double precision Wilson 
operator as two calls) of a) the relaxed SUMR inversion; b) the GMRESR{SUMR) inversion without deflation; c) 
the relaxed CG inversion, and d) the GMRESR(GG) inversion algorithm without deflation. 


steps while the eigenvectors are being calculated (to allow us to make use of a larger subspace). We 
deflated the eleven lowest eigenvectors of the Kernel operator to accelerate the computation of the 
matrix sign function (the Wilson operator tends to be well conditioned on a smeared dynamical 
overlap configuration, so this was all that was necessary). We tested the algorithm with 15, 30 and 
45 overlap eigenvectors. In all the inversions we used a Z 2 source vector, 6, where each component 
of the spinor was equally likely to be ±1. We used different choices of this source vector on our 
smaller test lattices, and it made no difference to the results. In all the plots for the inversions, 
we show the residual calculated during the CG or eigSUMR routine. For all the data, baring the 
issue mentioned in the discussion of [2] below, this does not differ from the true residual by any 
noticeable amount. For the GMRESR outer inversion, the true and calculated residuals did not 
differ significantly for any of our results. 

6.1. GMRESR and relaxation 

Figure [5] provides a comparison between the relaxed and GMRES preconditioned CG and 
SUMR routines, showing the gain of both the GMRES preconditioning and SUMR over CG. We 
generally will use the number of calls to the Wilson operator as a measure of the time for the 
inversion. Since the bulk of the cost of each algorithm is in the evaluation of the matrix sign 
function, and the cost of that is proportional to the number of times the Kernel operator is called, 
this is a reasonable measure to use, and it avoids various machine and run dependent fluctuations 
which would occur if we measured the total time. In all our figures, if the kernel operator was 
calculated in double precision, it was counted as two calls to the operator; while if it was calculated 
in single precision it was counted as one call to the operator. This compensates for the observation 
that a double precision Wilson operator is roughly half the speed of the single precision operator. 
Eigure [ 2 ] merely reproduces the results from earlier works. The curvature on the convergence of the 
relaxed GG and SUMR iterations is due to the relaxation. If we excluded the relaxation, the plot 
would be roughly a straight line, with the same gradient at seen on our plot in the early stages of 
the inversion. It can be seen that relaxation gives a vastly superior performance to the unrelaxed 
algorithm, and the GMRES algorithm (which we set to use five calls to the preconditioner) is 
about a factor of four superior to the relaxed algorithm. Thus, in total, these algorithms already 
give an order of magnitude improvement over the naive GG or SUMR algorithms. The small 
spikes seen in the GMRES(GG) curve are at the end of of the preconditioning steps. The plot 
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Figure 3: The convergence in terms of the number of Wilson operator calls (counting a double precision Wilson 
operator as two calls) of the a) undeflated GMRESR(relSUMR) routine; b) deflated GMRESR(deflated relSUMR) 
routine, with 30 overlap eigenvalues previously calculated to a high accuracy; c) undeflated GMRESR(relCG) 
routine, and d) the deflated GMRESR(deflated relCG) routine with 30 eigenvalues previously calculated. 


shows the measured residual; however, because of our aggressive relaxation strategy (we sacrificed 
a small amount of accuracy in the middle inversion for speed, since the GMRES inversion will 
correct for any small error), the true residual deviates from this by a small amount at the end of 
the inversion. The GMRES algorithm then resets the residual to its true value for the next call 
to the preconditioner. It can also be seen that SUMR is faster than GG for a single inversion, 
though by considerably less than the factor of two we might naively expect. This is partly caused 
by SUMR requiring a few more iterations to converge, and partly from the different relaxation 
strategies for GG and SUMR. 

6.2. Deflation 

We now turn to the effects of deflation. Eor these plots, we calculate and use thirty overlap 
eigenvectors during the deflation. We consider what occurs when we vary the number of eigenvec¬ 
tors later. We used five inversions to calculate the eigenvalues, which we label as calculating relCG 
or calculating relSUMR^ and then ran several inversions without improving the eigenvectors, and 
we label these runs as deflated relCG or deflated relSUMR. Eigure |3] compares the convergence 
of the deflated GMRESR(deflated relSUMR) and GMRESR(deflated relCG) algorithms with the 
undeflated algorithms. The plot is of the residual of the inversion against the number of calls to 
the kernel of the matrix sign function (given that almost all of the time required for the calculation 
is due to the computation of the matrix sign function, and the cost of the matrix sign function is 
proportional to the number of times the kernel matrix is called, this will be a good measurement 
of the overall cost of the calculation). Figure |4] shows the cost of the initial eigenvector calculation, 
plotting the number of kernel operator calls against the residual for those GMRESR(calculating 
relSUMR) routines used to calculate the eigenvectors, and comparing it against the undeflated and 
fully deflated GMRESR(relSUMR) routines. In figure [SJ we do the same thing for the iteration 
count of the inversion. Figures [6] and [7] repeat these results for the GG algorithm. There is 
some initial cost for the SUMR algorithm during the calculation routines before the eigenvalue 
computation. This is to improve the accuracy of , which we require to a reasonably accuracy 
while constructing the Galerkin projection otherwise the inversion algorithm will not converge to 
the correct result. This could have been achieved by keeping the accuracy of higher during 
the computation of the matrix sign function; however this will slow down the inversion algorithm 
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Figure 4: The convergence in terms of the number of Wilson operator calls (counting a double precision Wil¬ 
son operator as two calls) of a) the undeflated GMRESR(SUMR) routine; b) the deflated GMRESR(eigSUMR) 
routine, with 30 overlap eigenvalues previously calculated to a high accuracy; c) the first, third and fifth of the 
GMRESR(eigSUMR) routines where the eigenvectors where calculated. The long plateaus for calculating relSUMR 
1 indicates places where we needed to directly recalculate because of a poor estimate of 
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Figure 5: The convergence in terms of the number of calls to the matrix sign function of a) the undeflated GM- 
RESR(SUMR) routine; b) the deflated GMRESR(eigSUMR) routine, with 30 overlap eigenvalues previously calcu¬ 
lated; c) the first, third and fifth of the GMRESR(eigSUMR) routines where the eigenvectors where calculated. 
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Figure 6: The convergence in terms of the number of Wilson operator calls (counting a double precision Wilson 
operator as two calls) of the a) undeflated GMRESR(CG) routine; b) deflated GMRESR(eigCG) routine, with 30 
overlap eigenvalues previously calculated; c) the deflated GMRESR(eigCG) routine at various stages during the 
calculation of the eigenvalues. 



CG iterations 


Figure 7: The convergence in terms of the number of calls to the squared Dirac operator (containing two applications 
of the matrix sign function) of the a) undeflated GMRESR(GG) routine; b) deflated GMRESR(eigCG) routine, 
with 30 overlap eigenvalues previously calculated; c) the deflated GMRESR(eigCG) routine during the calculation 
of the eigenvalues. 
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considerably as more expensive matrix sign functions would be needed. We found that a better 
balance was achieved by keeping the accuracy of at the level needed to improve the eigenvectors, 
and have a small additional cost at the start of each inversion. 

On the configuration we used to generate these plots, the lowest eigenvectors of 75 I? had 
eigenvalues of ±0.114, while the twenty ninth and thirtieth eigenvalues were ±0.302. We therefore 
expect an approximately 0.302/0.114 = 2.65 gain in the iteration number for the inversion, if 
the deflation is successful. In practice, for the CG algorithm we see a gain of 265/98 = 2.70 in 
the iteration count (which is close enough to the optimum value), while SUMR gives a gain of 
329/165 = 1.99 (which is still far from the expected efficiency). The eigenvector residuals after 
five applications of the eigCG algorithm and used to deflate the CG inversions ranged from 0.01 to 
0.2, while the eigenvectors used in the deflated SUMR algorithm (also after five eigSUMR routines 
to calculate the eigenvalues) had residuals ranging from 0.0003 to 0.01. This tells us that firstly 
the eigSUMR algorithm is far more efficient at calculating the eigenvectors (which we expected, 
given the wrong Krylov subspace is used in the CG algorithm), and secondly that the deflation 
algorithm used in the eigGG routine works to its theoretical maximum even when the eigenvectors 
are exceptionally poor, while the projection/GMRES deflation used in SUMR did not work at full 
efficiently even with more accurate eigenvectors. We would require the eigenvectors to a greater 
accuracy, i.e. more inversions calculating the eigenvectors, to achieve the optimum improvement. 
Additionally, the full improvement of deflation is visible for CG after just the second eigenvalue 
calculation, while for SUMR the required iteration count gradually improves as the eigenvalues 
are computed. Nonetheless, despite the deflation method used for the CG algorithm’s greater 
efficiency, we still see that eigSUMR is more efficient than eigCG. 

We do not see quite the same gain when we consider the number of calls to the Wilson operator 
(in terms of the number of calls to the Kernel operator, the gain for GG was only a factor of 2.2). 
This is mostly due to the higher accuracy applications to the matrix sign function in the outer 
GMRES inversion. The deflation only affects the middle pre-conditioner, while the GMRES part 
of the algorithm remains untouched ~ a fixed cost regardless of whether we deflate. The deflation 
reduces the calls to the Wilson operator in the pre-conditioner at a similar rate to the reduction 
of the iteration count. The remaining inefficiency when converting iteration count into calls of the 
Wilson operator is due to the relaxation. 

Einally, we can see that the cost of calculating the eigenvectors is substantial, particularly 
for the SUMR routine where more calculating inversions are required before the gain from the 
deflation starts to become significant. The cost is usually largest for the first inversions, and 
even for the more mature calculations around a factor of 2-4 for GG and 4-5 for SUMR (i.e. one 
calculating inversion is twice to five times more expensive than a undeflated inversion). This is 
solely because we cannot relax the accuracy of the matrix sign function efficiently while calculating 
the eigenvectors. However, in a real-life simulation, only three or four of these steps are required 
(at least on the small lattices used in this study), while we will be able to use the deflated 
algorithm far more times. The gain from deflation will far outweigh the cost of the calculation of 
the eigenvectors. 

We also note that the deflation will become more effective as the quark mass is reduced, and on 
configurations with a non-zero topological charge (and thus exact zero eigenvectors of the massless 
operator). 

6.3. Eigenvector calculation 

We now consider the efficiency of the eigenvalue calculation, and whether our proposed relax¬ 
ation schema is successful. Eor this calculation, we used a standalone eigenvalue solver, starting 
from a Z 2 source, and continuing the eigenvalue calculation after the inversion had converged. We 
restarted the calculation after every twenty four diagonalisations of the lanczos vectors (i.e. when 
we were due to recalculate v^), or when the eigenvectors reached the required precision for that 
restart. At each restart, we performed an additional diagonalisation of the eigenvectors, using 
the computed eigenvectors and those from the previous two restarts. We used the final q vector 
from the previous run to restart the simulation, after re-orthogonalising it against the already 
calculated eigenvectors. This allowed us to further adjust the desired accuracy of the eigenvectors 
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as the calculation progressed, so we could always set the tolerance for the eigenvector run to be 
suitably lower than the expected residual at the next restart. We also restarted the eigenvector 
calculation when ||( 7 || — 1 grew larger than 0.1 (which indicates a breakdown in the algorithm due to 
the imprecision of the matrix sign function. If ||g|| grows too much then the algorithm can become 
unstable). This therefore does not fully replicate the situation in an eigSUMR inverter; however 
the rate of convergence is similar whether we use eigSUMR or this stand alone eigenvalue solver. 
We compare the residual for the eigenvector against the number of calls to the Kernel operator and 
the number of calls to the matrix sign function for both when we relax the accuracy of the solver 
according to the prescription described in section fS.ll and when we use a full accuracy matrix sign 
function for the whole calculation. In addition to the thirty eigenvectors, we used an additional 
80 Arnold! vectors when updating the eigenvectors. In an eigSUMR or eigCG calculation, this 
number would be restricted by the number of iterations needed for the inversion to converge. 

Figure [5] shows the residual of the eigenvectors plotted against the number of calls to the 
kernel operator for the first and last eigenvalues calculated, while figure [5] shows the residuals 
plotted against the number of recalculations of the eigenvectors. We see that, at least for the 
lowest eigenvector, the convergence is similar compared to the number of calls to the recalculation 
routine for both the relaxed and full accuracy eigenvalue solver. This means that our relaxation 
has not added to the number of iterations required to solve for the eigenvectors: the error in 
the eigenvectors is under control. When comparing the residue against the number of calls to 
the Kernel operator, we see a significant gain for the relaxed algorithm at the early stages of the 
calculation. 

However, while our eigenvalue routine works well to obtain the eigenvalues and eigenvectors to 
a low accuracy, it slows down considerably beyond a certain point, which depends on the number 
of eigenvectors calculated. If high accuracy eigenvectors are required in addition to the inversion, 
it may be advantageous to use eigSUMR to calculate the eigenvectors to a low precision, and 
then some other method such as inverse iteration or the Jacobi-Davidson algorithm to polish the 
eigenvectors to the required accuracy. The Jacobi-Davidson algorithm [d^, a more robust, faster 
and efficient generalisation of inverse iteration, generally converges well once the eigenvalues are 
known to a moderate accuracy, but less well if the eigenvalues are not known, although the cost of 
Jacobi-Davidson is also proportional to the number of eigenvectors required. The expensive part 
of the Jacobi-Davidson algorithm is a low accuracy inversion of the operator, which we already 
know how to do efficiently for overlap fermions. This is the opposite to what we have found for 
the Arnoldi process used in the eigSUMR routine. If high accuracy eigenvectors are required, we 
therefore suggest using eigSUMR to obtain an initial estimate of the eigenvectors and eigenvalues 
which can be used both as an input for the Jacobi Davividson routine and to construct a spectral 
preconditioner to accelerate it. 

We show the error \\v^ — Dv\\ plotted for the thirtieth overlap eigenvector in figure [TOl We 
see that the accuracy of remains slightly below the target accuracy, which suggests that our 
relaxation strategy works reasonably well. We calculate accurately after every twenty four 
diagonalisations of the Lanczos vectors, which causes the troughs seen in the plot. 

6.4- Varying the number of eigenvectors 

Finally, we consider the effect of varying the number of calculated eigenvectors, both on the 
convergence of the inversion and on the convergence of the eigenvector calculation. We consider 
15, 30 and 45 eigenvectors. We expect that the inversions will will improve as we deflate more 
eigenvectors, while it is less clear how the eigenvector calculation will fare. Figure [TT] shows the 
convergence of the deflated CG inversion as different numbers of eigenvectors are calculated, and 
Figure [12] shows the same thing for SUMR. Once again, the GG algorithm rapidly, after only a 
few iterations to calculate the eigenvectors, reaches the optimum convergence rate, while SUMR 
requires more than the five inversions we used to reach it. The ratio of the SUMR against GG 
improvement factor becomes worse as we deflate more eigenvalues, which means that the more 
eigenvalues we calculate, the more inversions are needed to have the eigenvalues calculated to a 
high enough accuracy for the deflation to take effect. For GG, this is not the case, and no matter 
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Figure 8: The convergence of the first and 30th eigenvectors in terms of the number of calls to the kernel operator, 
comparing the relaxed and full accuracy eigenvalue solvers. 
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Figure 9: The convergence of the first and 30th eigenvectors in terms of the number of calls to the matrix sign 
function, comparing the relaxed and full accuracy eigenvalue solvers. 
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Figure 10: The accuracy of when the accuracy of the matrix sign function was tuned by equation 1371 compared 
to the target accuracy given in equation 1331 

how many eigenvalues we calculate (at least up to the numbers we study here) only two or three 
inversions are required to calculate the eigenvectors to the required accuracy. 

Of course, the number of inversions required for the full improvement factor is one part of the 
cost of calculating the eigenvalues; the second is whether the cost changes per inversion for the 
CG and SUMR inversions used to calculate the eigenvectors. We display this in figures [T3] and 
m for the third calculating CG or calculating SUMR inversion, and in figures [15] and [H] for the 
initial calculating CG or calculating SUMR inversion. It can be seen that the additional cost per 
inversion required to invert extra eigenvalues is negligible. 

We conclude that in terms of computer time, the cost for the eigenvalue calculation for the CG 
algorithm is roughly independent of the number of eigenvalues calculated, at least when the size of 
the Krylov subspace used in each of the CG inversions is smaller than the number of eigenvalues 
needed. For SUMR, the computational cost to get the maximum gain for the deflation will increase 
as the number of eigenvalues increases. This difference is due to the improved preconditioning 
method used for the CG algorithm, which does not need the eigenvectors to be calculated to 
any significant accuracy to be effective; while the Galerkin projection method requires that the 
accuracy of the eigenvectors is comparable to the accuracy needed for the SUMR preconditioner 
in the GMRESR algorithm. With the number of eigenvalues we tested, deflated SUMR is still 
more efficient than deflated CG. That might change if more eigenvectors are included, or at lower 
masses, or fewer right hand sides are needed. Of course, for both CG and SUMR the required 
memory increases as the number of eigenvectors increases, and this may well be a limiting factor 
on some machines. 

In figures [T7] and [TH] we consider the convergence of the eigenvector calculation as more eigen¬ 
vectors are calculated. We consider the tenth eigenvector, which is calculated in all our runs, 
and the last eigenvector calculated in each run. It can be seen that the convergence of the tenth 
eigenvector improves as we increase the number of computed eigenvectors, while the final eigen¬ 
vector generally converges slower the more eigenvectors are calculated. The sharp drops in the 
eigenvector residual in figure [TH] occur when we reset the computation. They are partially caused 
by our diagonalisation of the eigenvectors against the eigenvalues computed from the previous 
iterations, which seems to particularly benefit the last few eigenvectors. If the eigenvectors are 
required in themselves, and not just for the deflation, then it may be advantageous to calculate 
more eigenvectors than necessary so that the convergence of the wanted eigenvectors improves. 
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Figure 11: The convergence of the CG inverter after deflating 0, 15, 30 and 45 eigenvectors. The maximum possible 
improvement over the original routine is given by the ratio of the lowest eigenvalue of 75 D not included in the 
deflation to the lowest eigenvalue of the operator. The ratio of the fifteenth to the first eigenvalue of 75 Z) is 2.1, of 
the thirtieth to the first eigenvalue 2.6 and of the 45th to the first eigenvalue is 3.1. The improvement in the number 
of iterations after deflation is a factor of 2.1 for 15 eigenvalues, 2.7 for 30 eigenvalues, and 3.1 for 45 eigenvalues. 
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Figure 12: The convergence of the SUMR inverter after deflating 0, 15, 30 and 45 eigenvectors. The improvement 
in the number of iterations for the deflated routines is 1.6 for 15 eigenvalues, 2.0 for 30 eigenvalues, and 2.2 for 45 
eigenvalues. 


28 













Figure 13: The convergence of the CG inverter while calculating 0, 15, 30 and 45 eigenvectors. 
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Figure 14: The convergence of the SUMR inverter while calculating 0, 15, 30 and 45 eigenvectors. 
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Figure 15: The convergence of the CG inverter while calculating 0, 15, 30 and 45 eigenvectors. 
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Figure 16: The convergence of the SUMR inverter while calculating 0, 15, 30 and 45 eigenvectors. 


32 














Calls to Kernel operator 



Iterations 


Figure 17: The convergence of the tenth eigenvector while calculating 15, 30 and 45 eigenvectors with the relaxed 
eigenvalue solver. 
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Figure 18: The convergence of the last calculated eigenvector while calculating 15, 30 and 45 eigenvectors with the 
relaxed eigenvalue solver. 


34 





















7. Conclusions 


We have tested the application of deflation methods to overlap fermions; focussing on the 
eigCG and eigSUMR algorithms. We construct a rigorous (though conservative) bound for the 
accuracy of the matrix sign function, and suggest a less rigorous and less conservative bound 
which seems to work well in practice. We see no loss in accuracy for the eigenvectors during the 
calculation. 

Our deflated inversion algorithms show the usual level of improvement. In particular, the 
spectral decomposition method we have used for the CG algorithm accelerates the inversions 
according to the theoretical maximum even for exceptionally low quality approximations of the 
eigenvectors. This pre-conditioner is impractical for other fermion formulations because of its cost, 
but for overlap fermions the cost is negligible. A method such as this should be preferred over 
the Galerkin projection for overlap inversions of D^D. For SUMR, we have to use a variation of 
the Galerkin projection, which performs less well until the eigenvectors are calculated to a high 
enough accuracy. 

As an eigenvector routine, our proposed eigGG algorithm performs particularly poorly, which 
is probably because the pre-conditioned inverter generates a different Krylov subspace than the 
optimal one for the calculation of the eigenvectors. However, it was still able to calculate the 
eigenvectors to a good enough accuracy to achieve the maximum possible acceleration from the 
deflation after only one or two right hand sides. The eigSUMR algorithm works well up to a certain 
accuracy, after which it shows a significant slowing down. To obtain high accuracy eigenvectors, 
it is beneficial to combine eigSUMR with some other routine such as inverse iteration or the 
Jacobi-Davidson method. 

There is a considerable cost in using eigSUMR or eigGG to calculate the eigenvectors over 
a straight inversion algorithm, because we are unable to relax the accuracy of the matrix sign 
function as aggressively as is possible in an inversion. This set-up cost will, however, be negligible 
compared to the gain if enough right hand sides are required. 

Because the deflation for the GG algorithm is more efficient than for SUMR, it remains an 
open question about which will be better in practical applications (if we do not also require high 
accuracy eigenvectors as well as the inverse - we expect that if the eigenvectors are calculated 
accurately enough then the Galerkin projection will work as well as the spectral preconditioner). 
Which of these methods is superior may vary from simulation to simulation. 

Obviously, deflation methods become less efficient on larger lattice volumes, since the increased 
density of small eigenvalues requires a larger number of eigenvectors to be calculated, leading to 
greater costs for the computational time and memory. However, given that our spectral pre¬ 
conditioner for the GG algorithm works with the maximum efficiency with even very low accuracy 
eigenvectors, and there does not seem to be an additional overhead for calculating more eigenvalues, 
we do not expect that there should be a significant increase in the computational cost at larger 
volumes (although, obviously, this needs to be confirmed). The memory requirement is likely to 
be more of a bottleneck, with the memory costs for the same degree of improvement likely to 
increase as O(U^), where V is the lattice volume, since the number of eigenvectors required for 
the same improvement is likely to increase with the volume. This should not be a problem on 
the relatively small (compared to other fermion actions) lattices which are currently accessible to 
overlap fermions, but it is something which will need to be addressed in the future. 
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Appendix A. Construction of the SUMR algorithm 

The Arnold! equation is UQ = QB, where i? is an n x n matrix, and Q is an A x n matrix 
containing the n orthonormal Arnold! vectors qi in its columns. This gives Q^U^UQ = B^Q^QB, 
Or fits = 1, i.e B is also unitary. This means that B can be decomposed in terms of Given’s 
matrices and a {7(1)” factor, which can be absorbed into the phase of the q vectors. Given that B is 
also upper Hessenburg (like all Arnold! matrices), the only Given’s matrices which can contribute 
to it are those with only diagonal and sub-diagonal terms. We can therefore decompose. 

Bn = Go(7o)Gi(7i) . . . G„_i(7n_i)G„(7n), (A-.l) 


where 
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, (A.2) 


where the f indicates the complex conjugate and am = a/I — |7mp, so that the Gm matrices are 
unitary. The final matrix is 


G„( 7 „) = diag(l, 1,... , 7 '„). 
By comparing the mth row, we can see that 

UQ 

m — Qm+1-^ m 1 


(A.3) 


(A.4) 


where Bm is the (m -I- 1) x (m) matrix defined in terms of the {m -f 1) x (m -f 1) Given’s matrices 
GiQ/i) and the (m -f 1) x m matrix G, 


We also know that 
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(A.5) 


(A.6) 


where (,) indicates the scalar product of two vectors. 

Our intention is to build up the Arnold! vectors qi via a short recurrence. We will first illustrate 
the process using an example when n = 4, and later generalise to arbitrary n. Equating equations 


36 



(lA.ll) and (|A.6I) gives 
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This gives, 

7i = -{qi,Uqi). 


The general form of this process is [34 1 


CTn =\/(l - \lu?) 
qn+l =Crnqn + llqn +1 


with qo = qo- We can construct the next q vector from equation (IA.4I1 . 

Uqm = (Qm+l)Go(7o)Gl(7l) ■ • • Cm- 
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With Qm+i = (go, 9i, 92 , 93 , ■ • ■), applying Qm+iGo gives 



Qm+iGo =(-7o9o +CTo9i,cro9o +7o'?l,'?2,93,---) 

=(9o, 91,92,93,94, ■ • ■)• 

(A.14) 

Similarly, 


Qm+lGoGiG 2 . . . G „_1 — ((9o, 9l, ■ • ■ , 9rt—l, 9n, 9™+!, Qn+2j ■ ■ ■) ■ 

(A.15) 

So we finish up 

with 



UQm = (■•■, 9m—1,9m, 9m+l)Gm- 

(A.16) 

We thus get 


U Qra — <^mQm+l 

(A.17) 

or 


Qm-\-l — {UQm “t” 

(A.18) 


This leads to the modified Arnold! algorithm shown in algorithm [TJ 


References 

[1] K. G. Wilson, Phys. Rev. DIO (1974) 2445. 

[2] P. H. Ginsparg, K. G. Wilson, A remnant of chiral symetry on the lattice, Phys. Rev. D25 
(1982) 2649. 

[3] P. Hasenfratz, V. Laliena, F. Niedermayer, The Index theorem in QCD with 

a finite cutoff, Phys. Lett. B427 (1998) 125-131. [arXiv:hep-lat/9801021], 

[doi:10.1016/S0370-2693(98)00315-3], 

[4] M. Liischer, Exact chiral symmetry on the lattice and the Ginsparg-Wilson relation, Phys. 
Lett. B428 (1998) 342-345. [arXiv;hep-lat/9802011] 

[51 H. B. Nielsen, M. Ninomiya, Absence of Neutrinos on a Lattice. 1. Proof by Homotopy Theory 
, Nucl. Phys B185 (1981) 20. 

[6] R. Narayanan, H. Neuberger, Chiral fermions on the lattice, Phys. Rev. Lett. 71 (1993) 
3251-3254. [arXiv:hep-lat/9308011], 

[7] H. Neuberger, Exactly massless quarks on the lattice, Phys. Lett. B417 (1998) 141-144. 
[arXiv;hep-lat/9707022] 

[8] H. Neuberger, Vector like gauge theories with almost massless fermions on the lattice, Phys. 
Rev. D57 (1998) 5417-5433. [arXiv:hep-lat/9710089], [doi:10.1103/PhysRevD.57.5417], 

[9] H. Neuberger, A practical implementation of the overlap-Dirac operator, Phys. Rev. Lett. 81 
(1998) 4060-4062. [arXiv:hep-lat/9806025] 

[10] P. Hernandez, K. Jansen, L. Lellouch, Finite size scaling of the quark condensate 

in quenched lattice QCD, Phys.Lett. B469 (1999) 198-204. [arXiv:hep-lat/9907022], 

[doi:10.1016/S0370-2693(99)01244-7], 

[11] J. van den Eshof, A. Frommer, T. Lippert, K. Schilling, H. A. van der Vorst, Numerical 
methods for the QCD overlap operator. 1: Sign- function and error bounds, Comput. Phys. 
Commun. 146 (2002) 203-224. [arXiv:hep-lat/0202025], 


38 




[12] A. Borici, A Schur Complement Approach to Chiral Fermions, PoS LAT2007 (2007) 065. 
[arXiv:0711.0508], 

[13] R. Narayanan, H. Neuberger, An Alternative to domain wall fermions, Phys.Rev. D62 (2000) 
074504. [arXiv:hep-lat/0005004], [doi:10.1103/PhysRevD.62.074504], 

[14] A. Borici, A. Kennedy, B. Pendleton, U. Wenger, The Overlap operator as a con¬ 
tinued fraction, Nucl.Phys.Proc.Suppl. 106 (2002) 757-759. [arXiv:hep-lat/0110070], 

[doi:10.1016/S0920-5632(01)01835-7], 

[15] R. G. Edwards, B. Joo, A. D. Kennedy, K. Orginos, U. Wenger, Comparison of chiral fermion 
methods, PoS LAT2005 (2006) 146. [arXiv:hep-lat/0510086], 

[16] A. Kennedy, Algorithms for dynamical fermions, in: proceedings of the ILFTN workshop 
‘Perspectives in Lattice QCD’, World Scientific, 2006. [arXiv:hep-lat/0607038] 

[17] Y. Shamir, Chiral fermions from lattice boundaries, Nucl. Phys. B406 (1993) 90-106. 
[arXiv:hep-lat/9303005] 

[18] D. B. Kaplan, A Method for simulating chiral fermions on the lattice, Phys. Lett. B288 (1992) 
342-347. [arXiv:hep-lat/9206013] 

[19] T.-W. Chiu, Optimal domain-wall fermions, Phys. Rev. Lett. 90 (2003) 071601. 
[arXiv:hep-lat /0209153] 

[20] N. Cundy, A. Kennedy, A. Schafer, A lattice Dirac operator for QCD with 

light dynamical quarks, Nucl. Phys. B845 (2011) 30-47. [arXiv:1010.5629], 

[doi:10.1016/j.nuclphysb.2010.11.017], 

[21] Z. Fodor, S. D. Katz, K. K. Szabo, Dynamical overlap fermions, results with hybrid Monte- 
Carlo algorithm, JHEP 08 (2004) 003. [arXiv:hep-lat/0311010], 

[22] T. A. DeCrand, S. Schaefer, Physics issues in simulations with dynamical overlap fermions, 
Phys. Rev. D71 (2005) 034507. [arXiv:hep-lat/0412005], 

[23] T. A. DeCrand, S. Schaefer, Chiral properties of two-flavor QCD in small volume 

and at large lattice spacing, Phys. Rev. D72 (2005) 054503. [arXiv:hep-lat/0506021], 

[doi:10.1103/PhysRevD.72.054503] 

[24] N. Cundy, et ah. Numerical methods for the QCD overlap operator. IV: Hybrid 

Monte Carlo, Comput. Phys. Commun. 180 (2009) 26-54. [arXiv:hep-lat/0502007], 

[doi:10.1016/j.cpc.2008.08.006], 

[25] N. Cundy, Current status of dynamical overlap project, Nucl. Phys. Proc. Suppl. 153 (2006) 
54-61. [arXiv:hep-lat/0511047], 

[26] N. Cundy, S. Krieg, T. Lippert, A. Schafer, Topological tunneling with Dynamical 

overlap fermions, Comput. Phys. Commun. 180 (2009) 201-208. [arXiv:0803.0294], 

[doi:10.1016/j.cpc.2008.09.010], 

[27] N. Cundy, S. Krieg, T. Lippert, A. Schafer, Dynamical overlap fermions with increased topo¬ 
logical tunnelling, PoS LAT2007 (2007) 030. [arXiv:hep-lat:0710.1705] 

[28] E. N. Zolotarev, Zap. Imp. Akad. Nauk St Petersburg 30 (1877) 5. 

[29] D. C. Sorensen, C. Yang, Accelerating The Lanczos Algorithm Via Polynomial Spectral Trans¬ 
formations, Tech, rep., Rice University (1997). 

[30] R. B. Lehoucq, D. C. Sorensen, Deflation Techniques for an Implicitly Restarted Arnold! 
Iteration, SIAM. J. Matrix Anal, and Appl. 17(4) (1996) 789-821. 


39 








[31] S. Hashimoto, et al., Lattice simulation of 2+1 flavors of overlap light quarks, PoS LAT2007 
(2007) 101. [arXiv:0710.2730] 

[32] G. Arnold, N. Gundy, J. van den Eshof, A. Frommer, S. Krieg, et al., Numerical methods for 
the QCD overlap operator. 11 Optimal Krylov subspace methods. Lecture Notes in Compu¬ 
tational Science and Engineering 47 (2003) 153-167. [arXiv:hep-lat/0311025] 

[33] N. Gundy, J. van den Eshof, A. Frommer, S. Krieg, T. Lippert, et al.. Numerical methods 
for the QCD overlap operator. 111. Nested iterations, Comput.Phys.Commun. 165 (2005) 

221-242. [arXiv:hep-lat/0405003], [doi:10.1016/j.cpc.2004.10.005], 

[34] Jagels, Reichel, A Fast Minimiml Residual Algorithm for Shifted Unitary Matrices, Numerical 
Linear Algebra with Applications 1(6) (1994) 555. 

[35] N. Gundy, S. Krieg, T. Lippert, Improving the dynamical overlap algorithm, PoS LAT2005 
(2006) 107. [arXiv:hep-lat/0511044] 

[36] A. Stathopoulos, K. Orginos, Computing and deflating eigenvalues while solving multiple 
right hand side linear systems in quantum chromodynamics, SIAM J.Sci.Comput. 32 (2010) 

439-462. [arXiv:0707.0131] 

[37] A. Kennedy, Fast evaluation of Zolotarev coefficients. Lecture Notes in Computational Science 
and Engineering 47 (2005) 169-189. [arXiv:hep-lat/0402038], 

[38] H. A. van der Vorst, C. Vuik, GMRESR: a family of nested GMRES methods, Numer. Linear 
Algebra Appl. 1 (4) (1994) 369-386. 

[39] 1. S. Duff, L. Giraud, J. Langou, E. Martin, Using spectral low rank preconditioners for large 
electromagnetic calculations, Int. J. Numer. Meth. Engng. 62 (2005) 416-434. 

[40] J. Erhel, F. Guyomarc’h, An Augmented Conjugate Gradient Method for Solving Consecutive Symmetric Positive I 
SIAM J. Matrix Anal. Appl. 21 (4) (2000) 1279-1299. [doi:10.1137/S0895479897330194], 

URL http;//dx.doi.org/10.1137/S0895479897330194 

[41] S. A. Kharchenko, A. Yu. Yeremin, Eigenvalue translation based preconditioners for the GMRES(k) method, 
Numerical Linear Algebra with Applications 2 (1) (1995) 51-77. [doi:10.1002/nla.1680020105], 

URL http;//dx.doi.org/10.1002/nla.1680020105 

[42] L. Giusti, C. Hoelbling, M. Luscher, H. Wittig, Numerical techniques for lattice QCD 
in the epsilon regime, Comput.Phys.Commun. 153 (2003) 31-51. [arXiv:hep-lat/0212012], 
[doi:10.1016/S0010-4655(02)00874-3], 

[43] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd Edition, Society for Industrial 
and Applied Mathematics, Philadelphia, PA, USA, 2003. 

[44] T. Chiarappa, K. Jansen, K.-I. Nagai, M. Papinutto, L. Scorzato, A. Shindler, C. Urbach, 

U. Wenger, I. Wetzorke, Iterative methods for overlap and twisted mass fermions. Computa¬ 
tional Science & Discovery 1 (1) (2008) 015001. [arXiv:hep-lat/0609023] 

[45] M. Luscher, P. Weisz, , Commun Math Phys 97 (1985) 59. 

[46] G. P. Lepage, P. B. Mackenzie, Viability of lattice perturbation theory, Phys. Rev. D 48 
(1993) 2250-2264. [doi:10.1103/PhysRevD.48.2250], 

URL http;//link.aps.org/doi/10.1103/PhysRevD.48.2250 

[47] G. Curd, P. Menotti, G. Paffuti, , Phys. Lett. B B130 (1983) 205. 

[48] G. L. G. Sleijpen, H. A. van der Vorst, A Jacobi-Davidson iteration method for linear eigen¬ 
value problems, SIAM J. Matrix Anal. Appl. 17 (1996) 401-425. 


40 






